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We propose a microscopic simulation for quark many-body 
system based on molecular dynamics. Using color confine- 
ment and one-gluon exchange potentials together with the 
meson exchange potentials between quarks, we construct nu- 
cleons and nuclear/quark matter. Statistical feature and the 
dynamical change between confinement and deconfinement 
phases are studied with this molecular dynamics simulation. 



At high baryon density, the nuclear matter is believed 
to undergo a phase transition to the quark matter be- 
cause of the color Debye screening and the asymptotic 
freedom in quantum chromodynamics (QCD) [Q. In 
qualitative estimates using the Bag model [D as well as 
the strong coupling lattice QCD Q predict a first order 
transition at baryon density (p) several times over the 
nuclear matter density (po — 0.17fm~^). However, re- 
alistic studies of the high density matter based on the 
first principle lattice QCD simulation are not available 
yet due to technical difficulties Q. In this situation, any 
alternative attempts are welcome to unravel the nature 
of high density matter. In particular, how the nuclear 
matter composed of nucleons (which are by themselves 
composite three-quark objects) dissolve into quark mat- 
ter is an interesting question to be studied. From the 
experimental and observational point of view, such tran- 
sition may occur in high-energy heavy ion collisions |^] 
and in the central core of neutron stars |^]. 

In this Letter, we propose a molecular dynamics (MD) 
simulation [^,^ of a system composed of many con- 
stituent quarks [||. As a first attempt, we carry out MD 
simulation for quarks with SU(3) color degrees of free- 
dom. Spin and flavor are fixed for simplicity, although 
there is no fundamental problem to include them. Time 
evolution of the spatial and color coordinates of quarks 
are governed by the color confining potential, the pertur- 
bative gluon-exchange potential and the meson-exchange 
potential. The confining potential favors the color neu- 
tral cluster (nucleon) at low density. However, as the 
baryon density increases, the system undergoes a transi- 
tion to the deconfined quark matter, since the nucleons 
start to overlap with each other. Our color MD sim- 
ulation (CMD) is a natural framework to treat such a 
percolation transition. The meson-exchange potential 
between quarks, which represent the non-perturbative 
QCD effects, helps to prevent the system to collapse. Al- 



though techniques are quite different, physical idea be- 
hind CMD with the meson-exchange potential is quite 
similar in spirit with the quark-meson coupling (QMC) 
model extensively used to study the nuclear matter from 
quarks pO| |. 

We start with a total wave function of the system as 
a direct product of single-particle quark wave-functions. 
The antisymmetrization is neglected at present. 
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Here A is the total baryon number of the system, (f>i is a 
Gaussian wave packet centered around with momen- 
tum Pi and a fixed width L. Xi is a coherent state in the 
color SU(3) space parametrized by four angles, ai,Pi,di 
and (pi. Although general SU(3) vector has six real pa- 
rameters, the normalization condition \xi\ = 1 and the 
unphysical global phase reduce the number of genuine 
parameters to four. Note that SU(2) spin coherent state 
parametrized by two angles has been used in the MD 
simulation of the many- nucleon system with spin ||ll|] . 

Time evolution of the system is given by solving the 
equations of motion for {Ri, Pi, at, Pi, 9i, tpi} obtained 
from the time-dependent variational principle 
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together with the classical Lagrangian, 
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dt 
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where H = {"S/ \ H \ '^). The explicit form of the equa- 
tions of motion reads : 
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As for the color-dependent quark-quark interaction, we 
employ the one-gluon exchange and the linear confining 
potentials. To take into account the essential part of 
the nuclear force, namely, the state independent short 
range repulsion and the medium range attraction, we in- 
clude the fj+Lo meson-exchange potential acting between 
quarks following ref. . The total Hamiltonian is writ- 
ten as 
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where — X"-/2 with A" being the Gell-Mann matri- 
ces, Vc is the confinement and one-gluon exchange terms, 
and Vm is the meson exchange term [Q. We introduce 
a smooth infrared cutoff to the confining potential in 
Vc{r) to prevent the long-range interaction beyond the 
size of the box in which we carry out MD simulations. 
We choose the cutoff scale Tcut = 3.0 fm, which is ap- 
proximately half of the length of the box. Typical val- 
ues of the parameters in the quark model for baryons 
read ||l2|, m = 350 MeV (the constituent-quark mass), 
as = 1.25 (the QCD fine structure constant), K = 0.75 
GeV/fm (the string tension). The meson-quark cou- 
pling constants ga{u))q are estimated from the meson- 
nucleon couplings ga(ui)N using the additive quark pic- 
ture: gaq = gcrAf/3 = 3.53 and g^q = g^jN/S = 5.85. 
The meson masses are taken to be /i^^ — 782 MeV and 
= 550 MeV. 

Some comments are in order here on the evaluation of 
the matrix elements H = {'^ \ H \ . 

(i) We have not taken into account the anti- 
symmetrization of quarks in the total wave function. Be- 
cause of this, the interaction between quarks in a color- 
singlet baryon is underestimated by factor 4 when one 
takes the matrix element of tftj. To correct this, we use 
effective couplings K''^ = AK and a'f^ = Aag throughout 
our CMD simulation. 

(ii) L (the size of the quark wave-packet) is chosen to 
be 0.35 fm (corresponding to the r.m.s. radius of the 
constituent quark of 0.43 fm). This is consistent with 
the typical value expected from the dynamical breaking 
of chiral symmetry [^3|. This value is to be used for 
taking the matrix element of the gluonic interaction Vc- 



On the other hand, the meson-quark coupling is intrin- 
sically non-local, since a and uj have their own quark 
structure. Besides, the meson-exchange interaction be- 
tween nucleons with the nucleon form-factor should be 
properly reproduced by the superposition of the meson- 
exchange interaction between quarks. To take into ac- 
count these facts, we use L®^ = 0.7 fm (corresponding to 
the r.m.s. radius of 0.86 fm) in taking the matrix clement 

ofVM- 

(iii) H = {'^ \ H \ '^) generally contains a kinetic energy 
originating from momentum variances of wave packets. 
However, when the width of the wave packet is fixed as a 
time-independent parameter, this kinetic energy is spu- 
rious and neglected in the present calculation. 

Let us now describe how to simulate the simplest three- 
quark system, namely the nucleon, in CMD. We first 
search for a three-quark state obeying the color neutrality 
condition 



(a = l. 
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This is satisfied by solving a cooling equation of mo- 
tion in the color space with a potential proportional 
to Y.i,o^i Ea=i(x»|A''|Xj)(Xi|A°|Xi) with random initial 
values of Xi- During this cooling procedure, the spatial 
coordinates of quarks are fixed, e.g. at the three corners 
of a triangle. 

If we start with three quarks in triangular position 
obtained above and kick each quark by a same amount 
of energy keeping the total momentum zero, the quarks 
start to have a breathing motion in 2-dimcnsional plane. 
Since the total color is conserved, the color-neutrality is 
maintained during this time evolution. 

By an initial kick to give the time-averaged kinetic 
energy of 74 MeV, the total energy of the nucleon become 
1269 MeV. Accordingly, the r.m.s. radius of the nucleon 
reads 0.46 fm in terms of L (which corresponds to the 
size of the quark-core of the nucleon) or 0.87 fm in terms 
of (which corresponds to the physical nucleon size 
for meson-exchange interaction). The "nucleon" here is 
certainly a semiclassical object which should be regarded 
as a mixture of the ground and excited states of three 
quarks. We use a collection of these nucleons as an initial 
condition for the CMD simulation of many quarks. Since 
the interaction among quarks in matter will eventually 
randomize the internal motion of quarks in the initial 
nucleon, the way how we kick the quarks does not matter 
for the final result. 

Now, let us study the phase change from the confined 
hadronic system to the deconfined quark matter. We 
simulate the infinite matter under the periodic boundary 
condition and see how the system responds to the change 
of the baryon density as well as to the energy deposition 
from outside. 
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To start with, nucleons constructed as above are ran- 
domly distributed in a box with the periodic boundary 
condition. At this stage, the total system is in its excited 
state. The minimum energy state of matter is obtained 
by the frictional cooling procedure, namely we solve a 
cooling equation of motion with frictional terms. During 
the cooling, spatial and color motion of quarks in the nu- 
cleon are artificially frozen, and the following equations 
are solved: 
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where hr and fip are damping coefficients with negative 
values and {i} means a set of three quarks in a nucleon 
to which i belongs. Under this cooHng procedure, the 
system approaches to a stable configuration with mini- 
mum energy. The system does not collapse due to the 
repulsive part of the meson exchange potential Vm- 

After the system reached its energy-minimum by the 
cooling, internal color and spatial motion of quarks are 
turned on and the normal equation of motion is solved for 
several tens of fm/c so that the system gets equilibrated. 
To study the excited state of the system, extra random 
motion is also given to the nucleons so that the system 
has a certain excitation energy. 

We judge the confinement/deconfinement by the fol- 
lowing criterion. If three quarks are within a certain 
distance dciustor and are white with an accuracy e, these 
quarks are said to be confined. This can be formulated 
as 
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All quarks are checked by this criterion without duplica- 
tions. The actual numbers we use are dciuster = 1 fm and 
£ = 0.05. 

Snapshots of matter in equilibrium for different exci- 
tation energies per quarks {E*) are displayed in Fig. 1. 
Quarks in the confined states are shown with thin col- 
ors and those in the deconfined state with thick colors. 
As E* increases, number of deconfined quarks increases 
as expected. However, some deconfined quarks still form 
three-quark clusters even for large E* . This implies that 
the deconfinement is caused not only by disintegration or 
percolation of clusters in the coordinate space but also 
by the color excitation inside each cluster. 
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FIG. 1. Perspective of matter with different excitation 
energies per quark E* at baryon density 0.78po • Particles with 
thin (thick) colors indicate quarks in the confined (deconfined) 
state. 

Figure 2 shows a "confined ratio of quarks", R = 
(number of confined quarks)/ (total quark number). At 
E* = 0, hadronic matter and the quark matter are well 
characterized by R although no sudden transition of R 
between the two phases is observed. For E* > 200 
MeV/q, R is less than 20% for all densities. 

To study the "thermal" property of the system, we fit 
the kinetic energy distribution of quarks by the classical 
Boltzmann distribution. Then, we can define an effective 
temperature T* for given E* . Note that T* is not really 
a physical temperature of the system, but is a measure 
of the averaged kinetic energy per quark. In Fig. 3, plot- 
ted is T* as a function of E*. For E* > 300 MeV/q, 
T* depends almost linearly on E* irrespective of baryon 
density. However, for E* = 100 - 200 MeV/q, T* for 
low-density matter increases rather slowly as a function 
of E* . In fact, this corresponds exactly to the region 
where the confined ratio of low-density matter changes 
in Fig. 2. This implies that, during the deconfinement 
process, the energy deposit from outside is consumed to 
melt the confined clusters (nucleons), which suppresses 
the effective temperature T* . 
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FIG. 2. Confined ratio of quarks as a function of E* for 
several different baryon densities. 



these reservations, the method proposed in this paper 
gives a starting point to study the statistical feature of 
the hadron-quark transition as well as to examine finite 
nuclei and the dynamics of heavy-ion collisions. Some 
preliminary simulation on the latter problem has been 
reported in 
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FIG. 3. Baryon density and excitation energy dependence 
of the effective temperature extracted from the Boltzmann fit 
of the kinetic energy distribution. 

In summary, we have proposed a color molecular dy- 
namics (CMD) simulation of the system with many con- 
stituent quarks. The system is approximated by the 
product of the wave packets with SU(3) color coherent 
state. Adopting the effective interaction between quarks, 
we study the transition from the nuclear matter to quark 
matter under the periodic boundary condition. At low 
baryon density (p) and low excitation energy (E*), the 
system is in the confined phase where most of the quarks 
are hidden inside the color singlet nucleons. However, 
as we increase p and/or E*, the partial deconfinement 
takes place due to the disintegration of color-singlet clus- 
ters both in the coordinate space and in the color space. 
This can be seen explicitly by the confined ratio and ef- 
fective temperature in Fig. 2 and Fig. 3. 

The results of this paper are still in the qualita- 
tive level. The refinement of interaction parameters, 
and the inclusion of flavor and spin degrees of free- 
dom as well as anti-quarks are necessary for more quan- 
titative discussions. The use of the antisymmetrized 
quark wave function is also an important future prob- 
lem p4| . The medium modification of the constituent- 
quark mass should be also considered in relation to the 
partial restoration of chiral symmetry. In spite of all 
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